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Hierarchy seems to pervade complexity in both hving and artificial systems. Despite its relevance, 
no general theory that captures all features of hierarchy and its origins has been proposed yet. 
Here we present a formal approach resulting from the convergence of theoretical morphology and 
network theory that allows constructing a 3D morphospace of hierarchies and hence comparing the 
hierarchical organization of ecological, cellular, technological and social networks. Embedded within 
large voids in the morphospace of all possible hierarchies, four major groups are identified. Two 
of them match the expected from random networks with similar connectivity, thus suggesting that 
non-adaptive factors are at work. Ecological and gene networks define the other two, indicating 
that their topological order is the result of functional constraints. These results are consistent with 
an exploration of the morphospace using m silico evolved networks. 
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I. INTRODUCTION 

Fifty years ago [P , Herbert Simon defined complex sys- 
tems as nested hierarchical networks of components or- 
ganized as interconnected modules. Hierarchy seems a 
pervasive feature of the organization of natural and ar- 
tificial systems [H The examples span from social 
interactions [U [5], urban growth [6l [7], and allometric 
scaling [5] to cell function [SHTTl [HI [T5] . development 
[13] , ecosystem flows [151 dS] > river networks [T7] , brain 
organization fTS' and macroevolution [THl 120] • But hier- 
archy is a polysemous word, involving order, levels, inclu- 
sion or control as possible descriptors |21| , none of which 
captures either its complexity or the problem of its mea- 
sure and origins. Although previous work using complex 
networks theory has quantitatively tackled the problem 
[3 m [^^H^ [?7H5T] , some questions remain: Is hierarchy 
a widespread feature of complex systems organization? 
What types of hierarchies do exist? Are hierarchies the 
result of selection pressures or, conversely, do they arise 
as a byproduct of structural constraints? 

A well established concept where such questions are ad- 
dressed involves the use of a morphospace |32H36j namely 
a phenotype space where a small set of quantitative traits 
can be defined as the axes. Here we take a step in this 
direction by combining morphospace and network theo- 
ries, taking the intuitive idea of hierarchy as the starting 
point: A pattern of relations where there is no ambigu- 
ity in who controls whom with a pyramidal structure in 
which the few control the many. Formally, the picture of 
hierarchy matches a tree of relations [37', ideally repre- 
sented by a directed graph. As shown in Fig. ([!]) the ele- 
ments of the system are represented by nodes connected 
by arrows establishing the map of relations of who af- 
fects whom. Accordingly, a measure of hierarchy should 
account for the deviations from this ideal tree picture. 
Such deviations occur because: (a) several elements are 



on the top, (b) downstream elements interact horizon- 
tally or (c) feedback loops are present. As shown in Fig. 
([T^), the tree-like picture matches the concept of genealo- 
gies, taxonomies, armies and corporations. Conversely, 
drainage networks in river basins |17| would define a re- 
verse anti-hierarchical situation, as depicted in Fig. ^p) . 
In this structure, multiple elements on the top merge 
downstream like an inverted tree. In between them, we 
can place a more or less symmetric web (Fig. [ij;) some- 
what combining both tendencies. But, in general, neither 
biological nor technological webs match the feedforward 
pattern. In most real systems, signal integration requires 
gathering inputs from different sources, while robust pro- 
cessing and control requires crosstalk and feedbacks, as 
represented by Fig. ([l]i), which are often organized in a 
modular fashion [55]. 

Hierarchy seems to pervade a coherent form of orga- 
nization that allows reducing the costs associated to re- 
liable information transmission |39j and to support effi- 
cient genetic and metabolic control in cellular networks 
[40j . However, a unified picture of hierarchy should not 
only provide a formal definition but also help understand- 
ing the forces that shape it. Here we provide the for- 
malization and the quantitative characterization of the 
morphospace of the possible hierarchies. The study of a 
large number of real networks and its comparison with 
model systems provide some unexpected answers to the 
previous questions. 



II. THE COORDINATES OF HIERARCHY 

Since hierarchy is about relations, our approach for- 
malizes the interaction between system's elements by 
means of a directed graph Q{V,E) [1] where Vi G V 
{i — is a node and {vi,Vj) £ E an arrow go- 

ing from Vi to Vj . This graph is transformed into the key 



2 



object for our methodology: the so-called node weighted 
condensed graph, Gc{Vc,Ec)- As Figs. ([IJi and [iji) 
show, Gc is a feedforward structure where cyclic mod- 
ules of G (the so-called strongly connected components 
{Si, Si}, hereafter SCCs) are represented by individ- 
ual nodes (thereby obtaining the condensed graph [1]). 
It is worth to mention that SCC detection has been 
shown as a powerful approach for subsystem identifica- 
tion [13 EH US] unraveling the presence of nested orga- 
nizations. In the case of a given Gc, every node has a 
weight ai that indicates the number of elements from G 
it represents (see section C 2 in the Appendix for details). 
Interestingly, such a graph represents the largest set of 
subgraphs of G (which, in this case, turn to be either 
isolated nodes or SCC^s) that can be properly ordered 
attending the causal flow defined by the arrows. This 
latter property will be the key to conceptually connect 
our definition of hierarchy with order. 

Our space of hierarchies will be a metric space fl, de- 
fined from three coordinates: treeness (T), feedforward- 
ness (F) and orderability (O), which properly quantify 
graph hierarchy. 



A. Treeness 

Treeness (T, being T in the range — 1 < T < 1) weights 
how pyramidal is the structure and how unambiguous is 
its chain of command. This measure covers the range 
from hierarchical (T > 0, Fig. ([l^)) to anti-hierarchical 
(T < 0, Fig. ^p)) graphs including those structures 
which do not exhibit any pyramidal behaviour [T = 0, 
Fig. ^)). As illustrated in Fig. ([l^-g), these graphs 
are characterized by taking the structure as a road map 
where we compare the diversity of choices we can take 
going top-down, i.e., following the arrows of the struc- 
ture, versus the uncertainty generated when reverting 
the paths going bottom-up. Such diversity is properly 
quantified using forward {Hf) and backward (Hh) en- 
tropies, respectively. Entropies are computed over di- 
rected acyclic graphs [21 13| . While section |d] of the Ap- 
pendix presents a rigorous description of these concepts, 
they can be briefiy presented departing from the node- 
weighted condensed graph GciVcEc) shown in Fig. . 
Notice that a directed acyclic graph is naturally a node- 
weighted graph and hence ensures G — Gc- From the 
graph Gcj we define two sets, namely M and fi. The first, 
composed by nodes with fci„ — 0, i.e. the set of maximal 
nodes and the second, the set of nodes with kout — 0, to 
be referred to as the set of minimal nodes. Now, let 11^;^ 
be the set of all paths starting in some maximal node. 
Since Gc is a graph without cycles, this set contains a 
finite number of elements tti, ...,tt[^. If Vi € M, the un- 
certainty associated to follow a given path starting from 
Vi and ending to some node in /i, h(vi) will be: 



where ¥{7rk\vi) is the probability that the path TTfe is fol- 
lowed, starting from node Vi € M. The average uncer- 
tainty we face to follow a path starting from some node 
in M win be Hf{Gc) = jM\J2v,eM H'^i) which is the 
general expression of the forward entropy of Gc- In a sim- 
ilar way but reversing the pathways a backwards entropy, 
Hb{Gc), of Gc can be obtained. Details of the derivation 
of Hb{Gc) are found in section [d| of the Appendix. 

Once the entropies accounting for the top-down and 
bottom-up path diversity generation are properly de- 
rived, we proceed to define a function with the aim of 
quantitatively grasp the deviations of Gc from the ideal 
tree picture of a hierarchical system. The explicit form 
is obtained by means of the normalized difference [^ of 
the two presented entropies, namely: 



Hf{Gc)-Hb{Gc) 
ma^{HfiGc),HbiGc)}' 



For the sake of consistency f{G) = if Hf{Gc) = 
Hb{Gc) — -occurring in the case where Gc is a lin- 
ear chain. The final value T{G) is computed as follows: 
First, let W{G) be set containing Gc and all subgraphs 
of Gc that can obtained by means of the application of 
a leaf removal algorithm (either top-down or bottom up, 
see section |d] of the Appendix for details). Then, T{G) is 
obtained by averaging / along all the members of yV(tJ): 



T{G) = if) 



w(e)- 



(1) 



In addition, T{G) = if Gc has no links -which can hap- 
pen, e.g., if G is totally cyclic. T{G) has been shown to be 
a very good indicator of the deviations of a given acyclic 
graph from the ideal hierarchical, tree-like picture [3]. 



B. Feedfordwardness 

In the Gc, since the elements within a SCC cannot 
be intrinsically ordered, SCCs constitute non orderable 
modules of the feedforward structure. As Fig. ([T|i) 
shows, SCCs represent a violation of the order at a given 
point of the feedforward condensed structure. Here, the 
size of the SCCs, but also their position in the feedfor- 
ward structure are key elements for the quantification of 
the impact of cyclic modules in the feedforward struc- 
ture, since the higher the SCC position, the larger the 
number of its downstream dependencies. According to 
this, we define feedforwardness {F, being < F < 1), a 
measure that weights the impact of cyclic modules on the 
feedforward structure of the graph, where cyclic modules 
closer to the top of G will introduce a larger penalty on 
hierarchical order than those placed at the bottom. For 
every path tt^ starting from the top of Gc we compute the 
fraction of the nodes that it contains against the actual 
nodes of G it represents. Formally, if f (tt^) is the set of 
nodes participating in the path tt^, 



h{vi) = - 



P(7rfe|v,)logP(7rfe|i;,), 



FiTTk) 
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regulatory networks ecological webs metabolism, language 



FIG. 1: Defining network hierarchy. The graphical representation of a tree-hke hierarchical graph (a), an inverted tree matching 
the anti-hierarchical graph (b), a non- hierarchical graph (c) with feedforward structures and a graph Q displaying cycles (d). 
In (e-f), all of them present 4 pathways (tti, 1x2,113, 7i"4) from maximals -M (top nodes)- to minimals -/i (bottom nodes)-. In (e), 
downstream diversity of paths is Hf{Q) = log 4, contrasting with the absence of uncertainty when reversing them Hi,{Q) = 0. In 
(f) the behavior is exactly the opposite: Hf{Q) = and Hh(G) ~ log4. In graph (g), the non-hierarchical feedforward structure 
presents the same forward and backward uncertainties, Hf{Q) = H{Q) = log 2. Chart (h) represents the node weighted 
condensed graph, Qc, computed by first detecting the SCCs and then collapsing every Si into single nodes. For 5*1 and ^2, 
the node weights are ai = 2 and 02 = 3, respectively. With a representative icon involving cones and balls (i-1) charts show 
TFO show the values for (a-d) graphs. In a vertical box, icons corresponding to three representative systems: gene regulatory 
networks (m), ecological networks (n) and metabolism and language (o). Icons are provided as an intuitive picture on how the 
fiow is organized within the network. They are used in fig. 2 for a more comprehensive view of the morphospace behavior. 
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where, as illustrated in Fig. ([T]), is the weight of node 
Vi in the node- weighted condensed graph. To obtain 
a statistical estimator of the impact of the location of 
cyclic modules within the causal flow described by the 
network, we first have to define the set Hm, which is the 
set containing all possible paths starting from the set of 
maximal nodes, M, and ending in any other node of Qc- 
Now, F{Q) is, thus, simply the average of F{-Kk) over all 
elements of IIm, i-c: 



F{Q) = {F)n, 



C. Orderability 



(2) 



As hierarchy is grounded on the concept of order, we 
need a descriptor that accounts for how orderable is the 
graph under study. Ranging from a fully cyclic graph to 
a feedforward structure orderability (O) lies in the range 
< O < 1 and it is defined as the fraction of the nodes 
of the graph Q that does not belong to any cycle. These 
nodes, therefore, make part of the fraction of the network 
that can be actually ordered. Such an orderable fraction 
provides a raw but very meaningful estimator about how 
ordered is the set of nodes within the graph (see Figs, 
([ij-l)). Formally, we define 0{Q) as: 



0(5) 



\{v^eVcnv}\ 

\v\ 



(3) 



Once we exposed the formal description of the differ- 
ent hierarchy indicators, i.e. treeness (T), feedforward- 
ness (F) and orderability (O), we proceed to collect them 
together to create a space fl where networks can be prop- 
erly evaluated and compared. 



III. THE DEFINITION OF THE 
MORPHOSPACE Q, 

Once our coordinates are defined, let us first see what is 
the repertoire of possible patterns that could be observed. 
According to our formalism, the hierarchical features of 
any directed network are given by a point u{0) in a 3D 
morphospace fl, being Q C [—1,1] x [0,1] x [0,1] (Fig. 
^) and Fig ([l7| of the Appendix). The point u{g) 
represents the graph Q by three coordinates, 



u{g)^{T{g),F{g),o{g)), 



(4) 



in the morphospace fi, according to its hierarchical prop- 
erties. Using the schematic representation of graphs 
outlined in Fig. ([T]), we define an intuitive icon asso- 
ciated to each kind of graph. As summarized in Fig. 
([2^), the perfect hierarchy is located at u{g) ^ (1, 1, 1), 
whereas the completely non-hierarchical system - a to- 
tally cyclic network- is located at u(C/) = (0,0,0). Inter- 
estingly, orderability and feedforwardness provide com- 
plementary information that defines forbidden regions 



of the morphospace. Since F{g) = 1 is only possible 
when 0{g) = 1, feedforward networks belong to the 
F{g) = 0{g) = 1 line. Given 0{g) ^ 1, not other 
F{g) 7^ 1 is permitted by definition. Attending to F and 
O, we find an interesting region, defined by 0{g) — 
and F{g) > 0. It is worth to stress that, while 0{g) 
and F{g) converge in their upper bound in a single value 
that defines the region of feedforward networks, they dif- 
fer when O goes to zero. This is because F{g) deals with 
the condensed graph while including condensed modules 
but 0{g) is about the nodes out of cycles. This little dif- 
ference allows to unravel a family of rare networks high- 
lighted by its extreme configuration -as shown by green 
icons in Fig. ([T^). Essentially they are formed by chains 
of small SCCs disposed in a feedforward structure. 

From these two axes accounting for the cyclic nature 
of networks, the coordinate T provides additional infor- 
mation about the organization of the resulting feedfor- 
ward structure after network condensation. Attending 
to the concept of a pyramidal structure, the plane sep- 
arating hierarchy from anti-hierarchy (T = 0) defines a 
family of symmetric structures, where the downstream 
path diversity is canceled by the uncertainty resulted 
when reversing these paths. Interestingly, when cycles 
are incorporated to this structure at T = 0, the result- 
ing structure is a bow-tie organization. In this particular 
structure, a large SCC occupies a central position in a 
quite symmetrized feedforward structure of inputs and 
outputs [45]. As shown in Fig. ([2^), the larger the SCC, 
the lower are the values of F and O. 



IV. NULL MODELS AND REAL NETWORK 
ANALYSIS 

Since null models do not consider optimal designs or 
functional constraints, they provide good insights about 
the part of the morphospace where no selection pressure 
is at work. Fig. ([2]d-c) shows that both homogeneous and 
broad random networks occupy (basically) the same re- 
gion in Q. This co-occupation, within the bow-tie plane 
with T{g) — 0, shows that random graphs appear lo- 
cated right in the middle between hierarchical and anti- 
hierarchical structures, independently of the type of de- 
gree distribution of the studied models. That implies, 
among other things, that they account for both bow-tie 
structures and feed-forward sparse webs. What about 
real nets? Here we use Nn = 125 real networks encom- 
passing 13 classes of natural and artificial systems (see 
SI2 for numerical details of u{g) values). As Fig. ([2ji) 
shows, a few isolated systems reach the boundaries of 11: 
a cell lineage located at u(C/) = (0.95, 1, 1) and a small 
social network located at u{g) — (—1,0.06,0.09). How- 
ever, most networks fall into four clusters. 

First, a group consisting of metabolic, neural, linguis- 
tic and some social networks is found at the lower part of 
the bow-tie domain, clearly embedded within the cloud 
of random graphs (Fig. ^) . Interestingly, random- 
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FIG. 2: The morphospace of possible hierarchies fl. (a) Different morphologies and their respective location within f2 (see 
Fig (1)). Green icons represent unlikely configurations, (b-c) The occupation of Q by an ensemble of random models. This 
set includes Erdos-Renyi (ER) graphs with different sizes (fOO, 250, 500) and average degrees (k) (see color bar). Symbols 
are proportional to network size, (c) Morphospace occupation of the Callaway growing network model overlapped with the 
ER ensemble as a reference. Three network sizes (100, 250, 500) and four connectivities {{k} = 2,4,6,8) are present (see also 
section E 1 of the Appendix for further models and other details), (d) The coordinates of the 125 real networks, colored and 
sized according to type and number of nodes, respectively. Cellular ([8] adapted from WormBase.org), Neuronal (^lOj taken 
from Newman's dataset collection). Social, Genealogy, Citations and Ownership were taken from Pajek database (Vladimir 
Batagelj and Andrej Mrvar http : //pajek. imfm.si/doku.php?id — data : index). Word corpora were generated according [16j . 
Food webs were taken Pajek dataset from Ulanowicz colecction. See references for Metabolisms |11H13| . Electronic circuits 
|16| . GRNs |14l 115) . Software TS' and Kinase networks \T^. Network data are available upon request. See section pi! 2 a| of the 
Appendix for more details about data compilation. Numbers and colors are network identifiers in SI2. CRN, MET, LANG, 
NEU, ECO and TECH stand for Gene Regulatory, Metabolic, Linguistic, Neuronal, Ecological and Technological networks. 
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ized nets of this group show a similar behavior although 
with a more central position within the cloud of random 
nets (see Fig. ([3|. An interesting case is given by the 
presence of bow-tie patterns in metabolic networks ^ 53j . 
They display a large central cycle, much larger than that 
observed in their randomized counterparts. This likely 
reflects the advantage of reusing and recycling molecules. 
The second group placed at the 0{G) ~ 1 plane shows 
a narrow band of feedforward nets including electronic 
circuits with —0.2 < T{Q) < 0.2 and software graphs 
slightly biased to negative T{Q) values. Here too the dis- 
persal seems consistent with what is expected from very 
diluted random graphs (see Fig. ([2j;)). This sparseness is 
a consequence of engineering practices focused on reduc- 
ing the wiring costs while keeping the system connected 

m- 

The third group displays slightly positive values of 
T{Q) and is composed of graphs with cycles of small 
size but with a predominant position at the feedforward 
structure giving rise to a very high O with variable F{Q). 
These are gene regulatory networks plus a protein kinase 
network. F{Q) dispersion is due to the variable size of 
modules located at the top of the structure. The spe- 
cial location in J7, far from the random cloud, is caused 
by a small fraction of genes, the DNA-binding elements 
(transcription factors) , located at the top of the network, 
which participate in cycles. Finally, the fourth group is 
defined by an isolated cluster of ecological flow graphs, 
located around u(^;) = (0.35,0.45,0.25). Their T{g) > 
values indicate a certain degree of pyramidal structure 
and the low 0(5)'s are consistent with an important role 
played by loops. The special status of these networks 
(not shared by other webs) is consistent with the well 
known picture of a trophic pyramid combined with the 
presence of recycling [111 IM] ■ 

The four clusters point to different scenarios pervading 
the origins of their hierarchical organization. One key ob- 
servation is that most data sets are found within the en- 
velope predicted by the ensembles of random graphs (see 
Fig. (|3|). Since these null models do not consider optimal 
designs nor functional traits, we conclude that, as it was 
reported in the context of modularity [SSI ISS] , hierarchi- 
cal order may be a byproduct of inevitable random fluctu- 
ations, which spontaneously generate graph correlations. 
Bow-tie networks, which have been suggested to define a 
flexible, robust and evolved type of systems [351 US j Per- 
haps resulting from selection for robustness would be a 
byproduct of the generation rules responsible for network 
growth. 



V. MORPHOSPACE ACCESSIBILITY BY 
DRIVEN EVOLUTION 

The previous results raise the question of how the voids 
in the morphospace need to be interpreted. In order 
to decide whether they are simply forbidden or have 
not been reached by evolution, we used an evolution- 




FIG. 3: Hierarchical order from fluctuations. Bow-tie net- 
works have been described as obvious examples of optimized 
structures that pervade robust functional behavior in cellular 
and technological networks. However, a close analysis reveals 
that the underlying organization of random networks with a 
given average number of links already displays such kind of 
pattern. Here a random network (a) with 200 nodes and an 
average degree of four links per node is built. The resulting 
Qc (b) displays a very clear bow-tie organization (c). The hu- 
man metabolic network (d), although non-homogeneous and 
resulting from evolution, displays a similar pattern (e-f). Both 
graphs are close within fl. Network layouts generated using 
Cytoscape software |58] . 



ary search algorithm |23| which, starting from a random 
configuration, tries to get the points of a evenly grid- 
ded partition of the morphospace. The results provide 
a picture about how accessible are the different regions 
of Q. Detailed information about the evolutionary al- 
gorithm is found in section |E 3| of the Appendix. Very 
briefiy, the evolutionary algorithm starts from a given 
set of small random graphs u e f2 belonging to the cloud 
of null models and networks. Given a target point u* of 
il, these graphs are evolved by a random process of link 
addition and deletion with selection of networks mini- 
mizing the distance ||uf — u*||. The number of nodes for 
every graph, \V\ = 50, remains constant in this process 
and graphs must belong to a single connected compo- 
nent. Only the number of links and their distribution 
are affected by the evolutionary algorithm. In this way 
the algorithm explores the network space, approaching 
the desired point and sometimes reaching it. The high 
computational cost of this experiment makes difficult to 
operate with larger networks. However, it is worth to 
note that their small size provides an advantage for the 
evolutionary search in the change of the network con- 
figuration since few changes in the connections have by 
general a deep impact in their structure. In this way the 
use of small network sizes contribute to an efficient ex- 
ploration, providing a coarse grained picture of network 
reachability. 

Fig. H reveals that the cloud of null models repre- 
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sented in Fig. ^p) is easily accessible, as indicated by 
the dark blue color of the region. As expected hierar- 
chical (T > 0) and anti-hierarchical regions (T < 0) are 
quite symmetric. Deviations are due to the dispersion 
produced by the finite-size statistics, specially when the 
resulting Qc become very little, by the imposition of a 
high number of cycles, as it happens for O = 0.15. This 
is the reason why, at O = 0.15, reachability is rather het- 
erogeneous. The solutions for this orderability are only 
possible by forcing to the network population to exhibit 
a large fraction of the nodes within cycles. This con- 
straint inevitably produces a Gc with just a handful of 
nodes. Then, as it happens for a broad number of topo- 
logical measures, values of T and F are extremely sensi- 
ble to small variations in network configurations. Such a 
trend is less dramatic when O increases (see = 0.5 and 
O — 0.85 of Fig. [4| since the fraction of nodes belonging 
to cycles is small enough to produce a rich combination of 
configurations in the resulting node-weighted condensed 
graphs. 

However, the most interesting results here concerns the 
presence of inaccessible regions, labeled in dark red color. 
Low levels of O seem to reduce the space of possible con- 
formations. At O = 0.15 the extremes of T and F are 
inaccessible under our in silico evolutionary experiments. 
Such a behavior is relaxed at O = 0.5 Here, a large re- 
gion of high reachability is observed for extreme values 
of T but, not occupied by real networks. This may indi- 
cate that a part of is accessible and yet not occupied, 
suggesting that the spontaneous correlations created by 
random fluctuations provide the source of order for free. 
As a consequence, non- adaptive processes would have 
played a major role in shaping hierarchies in nature |59j . 

Finally, an interesting trend is observed when O ap- 
proaches its upper bound. The larger O, the more re- 
duced is the range of the possible of F. Close to O = 1, 
there only exists graph configurations for F = 1, coincid- 
ing with feedforward networks encompassing, electronic 
circuits and software networks. 



VI. CONCLUDING REMARKS 

Both biological and cultural evolution operate under 
a number of deep constraints [BD| . Some of them result 
from the underlying rules of network growth and change, 
which strongly limit the repertoire of potential designs. 
An important question posed by evolutionary theory is 
the nature and relevance of such constraints in shaping 
the space of the possible. Our study provides a rationale 
for exploring the possible and the actual in complex net- 
works under a static view dominated by causal relations 
among components and modules. In this context, the 
inclusion of functionality, dynamics or weighted struc- 
tures has not being taken into account and should be the 
object for further work. By defining a general space of 
hierarchical webs, we are able to detect the presence of 
a rather limited domain occupied by real and random 




I 



■0.85 -0.65 -0.45 -0.25 -0.05 0.05 0.25 0.45 0.65 0.8! 




-0.85 -0.65 -0.45 -0.25 -0.05 0.01 



0.45 0.65 



FIG. 4: Morphospace accessibility by driven evolution of ran- 
dom networks. Here three sections of the morphospace were 



explored by an evolutionary algorithm (see section E 3 of the 



Appendix for algorithm definition) . To study the reachability 
of the regions of Q., the space was evenly spaced in 300 tar- 
get points u* distributed in a partition of 10 x 10 of the TF 
plane and three values of 0{Q) — {0.15,0.5,0.85}. Target 
points are labeled by crosses in the charts. Selection oper- 
ated during G — 10^ generations. Color indicates at which 
generation were acquired every grid in an average of 250 evo- 
lutionary experiments. Blue squares indicate that they are -in 
average- accessible after a few rounds of the process whereas 
red ones are unreachable before 1,000 iterations. 



null-models. The large voids surrounding these clusters 
of webs (defining four major groups) are partially inac- 
cessible and partially reachable, as shown by means of a 
directed evolution algorithm. The majority of webs dis- 
play a balance between integration of multiple signals and 
control over multiple targets under a bow-tie structural 
pattern. The computational nature of regulatory net- 
works and the combination of layers and cycles common 
to energy flows in food webs separate them from this large 
cluster. The matching of random and real webs in the 
first two clusters suggests that their hierarchical features 
can be accounted for from the spontaneous correlations 
associated to random graphs of a given degree, indicat- 
ing that the observed webs are simply the most proba- 
ble ones. By connecting network theory with theoretical 
morphology a powerful picture of complexity emerges, 
which allows us to both characterize hierarchical order 
and provide an evolutionary framework to explain how 
hierarchy emerges in nature. 
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Appendix A: Structure of the Appendix 

This Appendix presents in a self-contained way the 
conceptual issues that lead us to the rigorous formaliza- 
tion of hierarchy in complex networks. Although some 
of these concepts can be found in standard handbooks 
on graph theory -for example [1]-, we present all them 
from scratch in order to provide the reader with a con- 
sistent and clear mathematical apparatus. The aim is 
to remove any inconsistency in notation when introduc- 
ing original concepts. The type of object studied here is 
any kind of connected directed graph, i.e., any directed 
graph composed by a single component. Although the 
formalism could be applied over graph structures con- 
sisting in more than a single component -as usually hap- 
pens when working with random graphs- we assume that 
unconnected components have no causal relation of any 
kind and, therefore, the concept of hierarchy looses any 
meaning. As we shall see, the possibility to relate parts 
of the system by means of the topological information is 
basic for our definition of hierarchy. 

The appendix is structured as follows: First we present 
in an axiomatic way the features a given system based on 
relations among objects must hold in order to be consid- 
ered perfectly hierarchical. Deviations from this perfect 
configuration have to be properly quantified, giving rise 
to the three coordinates of the hierarchy. Then, we pro- 
vide the basis to work with, revising several basic defi- 
nitions of directed graphs and some of their most salient 
properties. We focus our attention on the condensation 
operation and on the layered structure of the resulting 
condensed graph. Basic definitions are followed by the 
proper definition of the hierarchy coordinates. Such def- 
initions are detailed and described with the aim of con- 
veying the reader the flavor and intuitions that underlie 
them. Once the formal framework is properly described 
we explore the regions of the defined morphospace occu- 
pied by model networks. This section is followed by the 
list of studied networks and the randomization methods 
we used to explore the relevance of the observed results. 
It is important to stress that the rigorous and systematic 
confrontation of real data with their randomized coun- 
terparts elucidates much of the possible origins of the 
observed patterns. Finally, we explore the space of pos- 
sible hierarchical configurations through an evolutionary 
algorithm whose rules are described in detail. 



Appendix B: The fundamentals of hierarchy 
(Postulates of Hierarchy) 

The perfectly hierarchical system will have the follow- 
ing three properties: 

• order, 

• reversibility and 

• pyramidal structure. 



Let us be more precise. Suppose we have a set ^ = 
{ai,...,a„} and a set relation R C A x A. Let T{R) 
be the transitive closure of R |27| , we will say that R is 
hierarchical if R has the following properties: 

• order, 

if {a,,ak) e T(R), then {ak,a,) i T{R) 

( There are no cyclic relationships: we can order 
the elements f28\j) 

• reversibility 

if {oi, flfe) G R and {oj, Ok) G R, then — aj 

( There is only one commander for any commanded: 
The chain of commands is reversible.) 



• pyramidal structure 
a) if {ai,ak) G 



R then, 



3a-i 



Hfc such that {ai,aj) ^ R{A commander commands 
more than a single element.) 

b) (3!afe e A) such that (Va^ e A\{ak}){{ak,ai) e 
T{R)) [There is only a single element which is not 
commanded by another element.) 

c) Let Ji be the set of commanders of element a;, 
i.e.: Ji = {ok & A : {ak,ai) e T{R)} If ai,aj e A 
are such that $ak G A such that {ai,ak) G i? or 
{aj,ak) G R, then \ Ji\ = \ Jj\. {All the elements of 
the bottom are subjected to a chain of commands 
of the same length. These last properties give the 
pyramidal structure of R) 

It is straightforward to realize that the directed, in- 
verted tree in which all the arrows go downwards starting 
from a single root naturally emerges as the graphical de- 
scription of the perfect hierarchical system -see fig. ([5|. 
The rigorous proof of this statement is given as a lemma 
in section |C 2| when the appropriate formalism is devel- 
oped. 




FIG. 5: A possible chart of relations depicting a perfectly 
hierarchical system. 



The measure of hierarchy must properly quan- 
tify the deviations from this perfect structure, 
considering these three postulated properties. 
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Appendix C: Graph definitions 

The whole set of relationships between elements of the 
systems under scrutiny define a directed graph. In any 
directed graph, links describe some kind of causal rela- 
tion between two elements of the set which is considered 
relevant for the system. We thus begin revising several 
facts about directed graphs. These basic notions can be 
found in any standard textbook on graph theory, like [T] . 



1. Directed Graphs 

Let Q{V,E) be a directed graph, being V = 
{«!,...,-(;„}, \V\ ~ n, the set of nodes, and E — 
{{vk^Vi), (vjjVi)} the set of arcs -where the order, 
{vk,Vi) implies that there is an arc in the following di- 
rection: Vk ^ Vi. Given a node Vi G V, the number 
of outgoing links, to be written as kout{vi), is called the 
out- degree of Ui, and the number of ingoing links of Vi is 
called the in-degree of Vi, written as kin{vi). The adja- 
cency matrix of a given graph t/, A.{Q) is a n x n matrix 
where Aij{Q) = 1 o {vi,Vj) e E; and Aij{Q) = oth- 
erwise. Through the adjacency matrix, kin and kout are 
computed as 

k^n{v^) = ^ Aj.iG); kout{vi) = ^^(^)- (^1) 



j<n 



3<n 



A path from node G y to node Vj G 7Th{vi,Vj) in 
a directed graph is an alternated sequence of nodes and 
links: 

TTk{Vi,Vj) = Vi, {Vi,Vj),Vj, ...,V£, {vi,Vj),Vj 

such that 

{v^,Vj), {vi,Vj) e E. 

(The k subscript is due to the possible presence of more 
than a single path from Vi to Vj). The length of a path, 
(■{T^k{vi,Vj)) is the number of edges appearing in the se- 
quence. We observe that some edges can appear twice or 
more in the sequence. In that case, we take into account 
all appearances of the edge. Since a path i^kivi, vj) is an 
alternating sequence of nodes and links, it seems natural 
to define two sets: v{TTk{vi,Vj)) as the set of all nodes 
present in the path and e{Trk{vi,Vj)), the set of all edges 
present in the path. Then we define the set 11(0) as the 
(possibly infinite) set of all paths that can be defined in 
G. 

A cycle in a directed graph is the subgraph formed by 
the edges and nodes defining a path which begins and 
ends in the same node, i.e, if 

Vi, {vi,Vj),Vj,...,ve, {vi,Ui),Vi 

such that {vi,Vj), {ve,Vi) E E. The set of cycles of 
a given graph is C(g) = {Ci,...,Ck}. If Ck G C(g) is 
maximal, i.e., 

$c, e c{g) : Ck c a, 



then, Ck is called a Strongly Connected Component, 
hereafter SCC. We will refer to the set of SCC's as 
Y,{Q) ~ {Si, Se}. Since every SCC is itself a sub- 
graph, we can refer to it as Sk = Sk{Vs^^, Es^), where 
J are the nodes and the edges of Sk , respectively. 
The underlying graph of a given directed graph 
Q, to be written as C/", is the undirected graph 
Q'"-{V,E'^) obtained by substituting all arcs of E, 
{vi,Vk) , {vj,Vs), ■■■■ by undirected edges giving the set 
E^ = {{vi,Vk},{vj,Vs}, ■■■}■ A directed graph Q is said 
to be connected if for any pair of nodes Vi,vi G V there 
is a finite, undirected path linking them; i.e., a finite se- 
quence 

Vi,{vi,Vk},Vk, .■.,Vs, {v„i,ve},ve 

where {vi,Vk}, {vm,vi} G Notice that now 

the links are depicted by unordered pairs, therefore 

{Vi,Vk} = {Vk,Vi}. 

Given a graph Q(V,E), a component is a maximal 
subgraph of Q by which, for every pair of nodes Vk,Vi 
belonging to it there is an undirected, finite path linking 
them. We refer to the set of components of Q as 

r(^) = {7i,...,7fc}. 

It turns out that a connected graph has only one com- 
ponent, the graph itself. In general, we will work with 
connected graphs. It is worth to note that, connected- 
ness may be lost during a randomization process but, as 
we shall see in section IE 21 the randomization methods 



we use respect the structure of components of the graph 
-i.e., connectedness if it be. 



2. Directed Acyclic Graphs 

A directed acyclic graph or feedforward graph - 
hereafter, DAG- is a directed graph characterized by the 
absence of cycles. The first consequence of the absence 
of cycles is the non-existence of infinite paths within -a 
finite- G- DAGs have been also referred to as ordered 
graphs because it is always possible to topologically sort 
a DAG [T]. A topological sorting consists in numbering 
all the nodes vi,V2, ■■■,Vn in such a way that arcs always 
point to nodes having higher numerical label than their 
origin, i.e., for any pair Vk,Vj such that k > j, then 

{vk,Vj) ^ E. 

This property will be crucial to justify the link between 
hierarchy, order and causal relation. This sorting im- 
plies that a DAG can represent a more or less entangled 
structure of interconnected causal processes. 

We stress in the analogy between order theory and 
DAG structures by defining following set: 



M = {u,&V: kUv,) = 0}, 



(C2) 
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to be named the set of maximal nodes of Q. Complemen- 
tarily, one can define the set of nodes ^ as 



{vi<^V : kout{vi) = 0}, 



(C3) 



which will be referred to as the set of minimal nodes of 

The set of all paths tti, tt^, from M to /i is indicated 
as IlMfiiG)- Given a node Vi E ii, the set of all paths 
from M to Vj is written as 



(C4) 



We observe that the length of the longest path of a finite 
DAG is always bounded. We will refer to this numerical 
value as L{g) which, taking into account the properties 
of the powers of the adjacency matrix, can be straight- 
forwardly obtained [29]: 

L{g) = max{fc : {3v,,v, G V : (A'=(g)),, ^ 0)}. (C5) 

A special type of DAG is to be defined now, since, as 
we will show below, depicted the perfect hierarchy. A 
directed tree in which all the leafs have the same length 
is a DAG with |M((7)| = 1 in which all elements but the 
one in M have kin = I7 either kout = or kout > 2 and 
in which \fvi,Vj G n{G), if Vk € M, then £{Tr{vk,Vi)) = 
£{Tr{vk, Vj j) -the leafs are all paths going from M to fi. 
See figure ([S]) for an example of this kind of graph. 

Let us now take the set-theoretic framework used at 
the beginning of this document and replace Ahy V (the 
set of nodes) and i? by £^ e V xV (the set of links) . Now 
we have a graph representation of the kind of relations 
that display a perfect hierarchical organization. It turns 
out that the type of graph defined above satisfies all the 
postulates defined in|B]is the graphical representation of 
a perfect hierarchical system. Let us state is as a formal 
lemma: 



Lemma: Under the hallmark described (see Postulates, 
section^B^, the kind of graph Q depicting the perfect hi- 
erarchy is a directed tree in which all the leafs have the 
same length. 

Proof: By the order property, Q has no cycles. By the 
reversibility property the graph nodes of the graph have 
either kin = or kin = 1- By the pyramidal property a) 
the nodes have either k^ut — or k^ut ^ 2. By pyramidal 
property b) the graph is connected and |M(CJ)| — 1 and 
by pyramidal property c), if Vk G M, then Vw^, Vj S fJ-iQ), 
then £{T:(vi^,Vi)) = £(Tr{vf;,Vj)) (notice that there is only 
one path to go from Vk G M{g) to a given Vi e fJ.{G))- 
Therefore, Q is directed tree in which all the leafs have 
the same length. 

Now that we have the mathematical characterization 
of the perfectly hierarchical graph, we restart the concep- 
tual presentation that will help us to properly evaluate 
deviations from the ideal graph. 

a. Condensation: Obtaining a DAG from any directed 
graph 

Now we consider again the wide class of directed 
graphs, not only DAGs. Let T,{Q) = {Si, ...,Sk} be the 
set of SCGs. We will build the condensed graph of G, 
to be referred to as Qci in which every SCG is merged 
into a single node, maintaining the links that connect a 
node of such SGC with nodes out of such SGC. More 
formally. 



Vr 



s(a)U 




i.e., nodes of Qc are either SGCs of Q or nodes that do 
not belong to any SGC of Q. Consistently, edges on Eq 
connect nodes of Vq in such a way that 



{v,,v.j) e E 

(yv„v, e Vc){v„v.j) eEc^ { {v^ = Ske ^{g),v^ eV)A {3ve G St) : {{vi,Vj) G E) 

{v, = 5fe G ^{Q),vj ^Si& E(5)) A {3vm G Sk,3vs G St) : {{v„„Vs) G E); for k ^ £. 



It turns out that, by definition, Qq is a DAG. 



b. Node-weighted condensed graph 



Below we have an example of the condensation oper- 
ation, step by step. We have a graph Q (left). We then 
identify the SCGs of Q (center) and collapse all nodes 
belonging to a SGC in a single node, thus obtaining Qc 
(right): 



A conceptual step beyond the condensed graph is the 
node-weighted condensed graph. All the computations 
presented in this paper will be performed over such a 
graph. Let us formally describe it: In this graph, a weight 
ai is assigned to every node Vi G Vc in the following way: 



Notice, as illustrated in figure ^ that Qq is a DAG, 
as expected. 



\S,\^v, G S(g). 
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V 



FIG. 6: The three stages of the condensation algorithm. 

Therefore, the sequence of 

A = ai, ■■.,a\vc\ 
will be the sequence of weights of Qc ■ Clearly, 



*<|Vc| 



^\V\. 



In plain words, ai accounts for the number of nodes be- 
longing to V represented by a single node in V^. In the 
forthcoming example we detail the obtaining of such a 
graph from a given directed graph. We have the graph 
-left-; then, we identify its SCCs -center- and, then, - 
right- we label the nodes of Qc with their corresponding 



a s: 




FIG. 7: Obtaining the node-weighted condensed graph Q. 



c. Dissection of the layers of the graph 

The objective of this section is to rigorously define the 
dissection in layers of a given DAG. The non-cyclic na- 
ture of such graph structures enables us to define a fi- 
nite set of layers. Layers can be identified through a 
backward or bottom up leaf removal (LRf,) algorithm or 
through a forward or top down leaf removal (LRf) algo- 
rithm [30]. Although LRf and LRi, algorithms identify 
the same number of layers, the specific composition of ev- 
ery layer is not necessarily the same. In addition, using 
these two algorithms, we can build two sets of subgraphs 
by successively removing layers of the DAG. As we shall 
see, these two sets will be very useful to derive the pro- 
posed measure of hierarchy. 

Accordingly, we define two partitions of Vc, W = 
{wi , . . . , uJm } and W = {wi , . . . , ai m}- The number of layers 



relates to L{Qc) -see equation (C5|- as follows 



\w\ = \w\^L{gc) + i. 

The members of such partitions are the nodes defining 
the layers of the DAG, computed by either a or 
LRf algorithm -depending on which partition we gener- 
ate, either W or W , respectively. Specifically, using the 



bottom up approach -i.e., using a LR^ algorithm- the 
first member of such partition is defined as the following 
subset of nodes: 

cji = {vi e Vc : Kut{vi) = 0} 

and, using the top-down approach -i.e., using a LRf 
algorithm- the first member is defined as, 

uJi = {v, G Vc : kin{v,) = 0}. 

Clearly, cji = /i and uji = M. With the above subsets of 
Vc we can define the graphs Qi{Vi,Ei), and Gi{Vi,Ei) 
in the following way: 



Vi = Vc \ uji; El = Ec\ {{vi,Vk) : Vk G uji} 



and 



Vi^Vc\ wi; Ei=Ec\ {{v^, Vk) : Vi G cDi}. 

respectively. Similarly, we build W2 and 0)2 as: 

W2 = {vi G Vi : kout{vi) = 0}, 
W2 {vi G Vi : hn{vi) = 0}. 

which, in turn, enables us to derive V2, i?2, V2 and E2: 

V2 = Vi\ UJ2; E2= Ei \ {{vi, Vk) : Vk G UJ2} 
V2^Vi\ UJ2; E2^ Ei\ {{v„ Vk) : Vi G (iz}- 

In the general case, if Vq = Vq = 

Vt = Vi-i \ LOf, Ei = Ee^i \ {{v^,Vk) : Vk G uJe}- 
Ve = Vi-i \ Qi; Ei = E^^i \ {{vi, Vk) : Vi G cj^}, 

where 

i<i<Ligc). 

It is worth to note that 
and that 

The two previous sequences of subgraphs can be ordered 
by inclusion, namely 



and 



^i(ec) c ••• c cJi c cJc, 



SLiGc) c ... c Gi cGc- 



The dissection of Gc in the above described two sequences 
of subgraphs will enable us to exhaustively explore the 
role of all layers in the further hierarchy measure. We 
finally define the set yV{G), containing the graph Gc and 
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all non-empty subgraphs obtained by means of the appli- 
cation of a leaf removal algorithm (either bottom up or 
top down) and which contain at least one link: 

W{g) - {gcGLigc)-!^ SuGLiGc)-!. Gi}. (C6) 
It is not difficult to check that 

\w{g)\^2L{gc)-i. 

Let us provide an example in figure Q: Starting from 
a DAG (left), we identify the layers uji,...,uj5 using a 
LRb algorithm (center), and the layers (2;i, cl^s using 
a LRf algorithm (right): Furthermore, we can obtain 





FIG. 8; The LR algorithms and the identification of layers. 
In this case, IVKj = 5 and L{g) = 4. 



the sequence of graphs ordered by inclusion G^Gi, ■■■^G'k 
-notice that the last graph consists of two isolated nodes: 



FIG. 9: The sequence of subgraphs G ---jGi- 



And the sequence of graphs ordered by inclusion 
G, Gi, Gi- 




FIG. 10: The sequence of subgraphs QyGi, ■■.,^4: 



Appendix D: Detailed derivation of the Coordinates 
of Hierarchy 

This section is devoted to the formalization of a quan- 
titative descriptor of hierarchy. Such a descriptor is de- 
fined from three hierarchy coordinates. These three coor- 
dinates are: Orderahility (O), Feedforwardness (F) and 



Treeness (T). Thus, any directed graph will define a 
point in such space, described as follows: 

u(g)^(T(g),F((?),0(g)). 

From the coordinates of u(Q) we can extract information 
about the structure of the net under the conceptual back- 
ground of hierarchy presented in section [B] based on the 
concepts of order, reversibility and pyramidal structure. 
Through the values of these coordinates we must be able 
to identify and properly quantify deviations from such an 
ideal behavior. As we shall see, such three components 
naturally arise as long as we go in depth in our inquiry for 
a hierarchy estimator. The order x = T,y = F,z = de- 
rives from clarity issues related to the visualization. For 
the sake of clarity in the exposition, however, the order 
of the sections in which we present such measures will 
not follow such T, F, O structure. Instead, we define the 
coordinates in such a way that we go from the simplest 
(O) to the most complex one (T). 

We are now ready to define the coordinates of hierar- 
chy. 



1. Orderahility, O 

This is the z coordinate ofu{Q). 

The Orderahility, O, of the graph Q is the fraction of 
nodes not belonging to any cycle -which are, by defini- 
tion, non-orderable structures. In a more formal way, 
let C/ be a directed graph and Gc be its condensed coun- 
terpart. The orderahility of the graph 0{G) is defined 



OiG) 



\{v^evcnv}\ 

\v\ 



(Dl) 



In terms of the node-weighted condensed graph, we can 
rewrite the above expression as: 



OiG) 



\{v, e Vc : g., = l}| 
1^1 



Let us provide an example: In the graph depicted be- 
low (left), \V = 7|, \Vc\ = 5 (right). The grey circle 
depicts a condensed SCO. \{vi £ Vc H V}\ — 4, there- 
fore, 0{G) = 4/7: 




FIG. 11: Computing the orderahility of the graph. 

This is a raw indicator of the fraction of the net that 
can be ordered. A tree or, more generally, a feed-forward 
network will display 0{G) = 1, in agreement to their 
orderable nature. In contrast, a network which is totally 
cyclical would display 0{G) = 0. 
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2. Feedforwardness, F 

This is the y coordinate ofu(Q). 

Beyond the nodes than can be ordered, now we want 
to know the impact of the non-orderable regions of the 
graph over the potential causal paths described by it. In 
raw words, where, within the causal flow, we find the 
non-orderable regions. This is captured by the to-called 
FeedForwardness, F[Q), a measure centered on the paths 
of the graph Qc beginning in the set of maximal nodes, 
M -see equations (C2 C3| and (C4|. Specifically, for 
every path going from M to Vc\ M we define a 
function, F, which evaluates the quotient between the 
number of nodes and the overall weight of the path. To 
cover all nodes of the graph, the numerical value of the 
coordinate will be averaged over the sequence of graphs 
obtained when applying a algorithm. If Qc is finite 
both M and V\M are finite and we can safely compute 
averages of these observables. 

To put the things in a more concrete way, let Q be 
a graph and Qc its condensed counterpart. If TTfc e 
^Mfi{Qc)i we define the function F of ttj,, F[-nk), as: 



l«(i'fc)l 



For example, if we take the graph studied above and we 
choose the highlighted path, to be named tt^ : 



we have that F{TTj) — |. 

To compute the overall average we then build the func- 
tion g: 

g{Q) ^ ^('^^)' 



and we compute 5(^1), 5'(^L(e)-i) -see section C2c 
for the definition of ^^'s: 

9{Qk) = E ^('^») 

We finally average over all CJ's and Qc, skipping QL{g) 
because, by its very definition, this graph has no paths: 

F{Q) = 9{Gc) + Ek<Ligc)aiSk) ^^^^ 
\^MtJ.iQc)\ + Efc<L(ec) \^Mt,{Qk)\ 

And, for the sake of consistency, if the graph consists in 
a single node. 



11^1 



1) 



(FiQ) = 0). 



The interested reader can check that, with equation (D2 1, 
we cover all paths from M to y \ M. 

Again, we will find F{Q) = 1 in DAGs, as well as 
F{Q) = in networks consisting of a single strongly con- 
nected component. The combination of 0{Q) and F{Q) 
tells us how are cycles located within the net and the 
impact they have in the order of the structure. The com- 
bination of 0{Q) and F{Q), thus, provides us interesting 
information on how the net is globally organized. 



3. Treeness, T 

This is the x coordinate ofu{Q). 

This measure accounts for the reversibility and the 
pyramidal structure. The forthcoming coordinate is of 
information-theoretic nature and can only be computed 
over the condensed graph Qc , since it can only be defined 
over DAGs. The impact of cycles (non-orderable regions 
of the graph) is already evaluated by the two previous in- 
dicators. This new coordinate largely relies on the com- 
putation of a special kind of entropies which, for the sake 
of brevity, are not going to be derived in detail here. We 
refer the interested reader to [21 [3] . 

Intuitively, we first observe that deviations from re- 
versibility can be properly captured measuring the uncer- 
tainty in reversing a given path, a measure called topo- 
logical reveresibility [3]. This measure of uncertainty is 
provided by the statistical entropy over the set of paths 
present in a given graph. It turns out that deviations of 
the pyramidal structure can be also quantified through 
an information-theoretic measure, namely, from the dif- 
ference between the amount of statistical entropy gen- 
erated when we cross the graph according to the flow 
direction depicted by the arrows and the uncertainty in 
reversing the paths [3]. The spirit of the measure is to 
compare the creation of alternatives -new paths- in a top 
down exploration of the DAG against the irreversibility 
of these paths. If the quantification of the uncertainty in 
reversing the paths displays a lower value than the quan- 
tification of alternatives, we say that, qualitatively, the 
net shows a pyramidal structure. If the creation of new 
paths in a top down exploration of the graph is quanti- 
tatively equal to the uncertainty in reversing them, then, 
there is no argument to attribute to the net any pyra- 
midal shape. Finally if the situation is opposite to the 
former one, the graph will display a funnel-like or in- 
verted pyramid structure. Below we present the above 
mentioned entropies in more detail. 



a. Backwards entropy 

Let Q he a, directed graph and Qc its condensed coun- 
terpart. Let A{Qc) be the adjacency matrix of Qc. We 
first define the \Vc \ M\ x \Vc \ M\ matrix B{Q) in the 
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FIG. 12: We have a graph (left). Imagine that we know 
that some causal chain ended in node a. We have 5 different 
paths (center, inside the rectangle) to reach a from the top 
-the maximal nodes. So there is uncertainty in reversing the 
path. Finally (right), we observe that there is no uncertainty 
in reversing the causal path that ended in /3, since there is 
only one path to go to /3 from the maximal nodes. In this 
graph Hb > due to the uncertainty arisen when reversing 
the paths that ended in a. 



FIG. 13: We have a graph (left). Imagine that we know that 
some causal chain begun in node /3. We have 4 different paths 
(center, inside the rectangle) to reach a node belonging to the 
set of minimal nodes. So there is uncertainty in following the 
path. Finally (right), we observe that there is no uncertainty 
in following the causal path that begun in a, since there is 
only one path to go from a to the minimal nodes. In this 
graph, Hf > 0, due to the uncertainty arisen in following the 
paths starting in /3. 



following way: 



iyv,,Vj : = 1) B{g), 



(D3) 



and {yvi,Vj : Aij = 0) B{Q)ij — 0. From this defini- 
tion, we obtain the explicit dependency of the probabiUty 
of crossing Vj when we start to reverse a path from w^, 
V{vj Vi). Thanks to the Z? AG-like nature of the con- 
densed graph, we can compute such a probability directly 
from the powers of the adjacency matrix [2], namely, 



-0= E 

i<k<L{gc) 



(D4) 



Now we compute the average amount of uncertainty we 
have to face when reversing a path. Specifically, we have 
path starting at some node in M and ending at node 
Vj G fi. We want to know the uncertainty of recovering 
this path if we go backwards, i.e., from the node Vj € fito 
a given node in M. This average uncertainty is provided 
by the following entropic functional [3]: 

Hb{gc) = ^^J2 J2 P(i'. 4-Wfc)-logfc„(z;fe). (D5) 



b. Forward entropy 

Now we compute the forward version of the above de- 
fined entropy. In this case, we need to compute the prob- 
ability to cross node Vk departing from Vi G M according 
to the causal flow -not reversing it, as above. The ex- 
plicit expression of this probability is defined from matrix 
B'iG): 



{\fv„Vj -.A,, = 1) B'iG), 



and {Wvi^Vj : Aij = 1), B'{Q)ij — 0. Then, in analogy to 
what we found above -equation (D4|-, we have that the 



probability to cross node Vj if we started a path in Vi is 
nv.^v,)= E (D6) 



i<k<L[gc) 



The average of the amount of uncertainty emerging when 
we want to follow a path from which we know that starts 
in a given Vi e M and ends in a given node of /i will now 
be ij: 



(D7) 



where M is the set of maximal nodes 



c. Treeness 




FIG. 14: Gaining intuition: Why do Hb and Hf grasp the 
hierarchical nature of graphs. The graph on the left displays 
Hf — Hb > and it is thereby hierarchical. On the contrary, 
the graph on the center displays Hf — Hb < 0, being thus 
antihierarchical. Finally, on the right we have a graph such 
that Hf — Hb ^ 0, thus it is neither hierarchical nor anitihier- 
archical. 



The presented coordinate compares the creation of new 
paths following the causal flow defined by the arrows 
against the uncertainty to reverse them. The creation 
of new paths is due to the existence of more than a sin- 
gle alternative to leave a given node. If we create new 
paths we are thus creating information -quantified by Hf. 
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However, such information can be destroyed by the un- 
certainty in reversing the paths -evaluated by Hh. In- 
tuitively, a graph having pyramidal structure will have 
the balance Hf — Hi, positive, whereas a graph having 
inverted pyramidal structure will display Hf ~ H^ nega- 
tive. A completely random DAG displays H f — H\j k, Q 
[3]. It turns out that, properly manipulated, such infor- 
mation measures are the perfect indicators of the pyra- 
midal properties of the graph, a key ingredient of our 
hierarchy coordinates. 

In order to generate a normalized estimator (between 
— 1 and 1) accounting for the balance between Hf{Qc) 
and Hi,{Qc) we define f{Q) as follows: 



, . ^ Hj{gc)^H,[gc) 

^^^> max{Hf{gc),Ht{gc)}' 



(D8) 



The treeness coordinate of a DAG, to be indicated as 
T{g), will be the average among the i(tJc) ^ 1 sub- 
graphs tJi, ...,t/L(ec)-i> the L{gc) - 1 subgraphs 
gi,...,gfc,---,gLfgr-|-i and g Itself, i.e., the set W{g) de- 
fined in equation (C6| [32l, i.e.: 



ng) 



1 



2L{gc) - 1 



fiG) 



f{g^) + ^{g^ 

KLiGc) 



I We I 



(D9) 



e.6W(e) 

= {.f)Wg, 

where Wg is the set of subgraphs of C/c obtained through 
the application of a lead removal algorithm (either bot- 
tom up or top down) and which contain at least one link, 
as defined in equation (C6). 



FIG. 15: The need for detailed layer analysis. All four graphs 
display Hf > and Hf, — 0. However, the two graphs on 
the right contains a violation on the hierarchy assumptions, 
namely, the need for pyramidal structure at all levels. The 
dissection allows to identify and penalize such deviations. 

For the sake of consistency we explicitly define the 
behavior of some limit, potentially problematic cases. 
When obtaining the sequence of g's and ^'s the graph can 
break into more than a single connected component. Let 
r(CJi) ~ {7i,...,7fc} the set of components of our graph 
C/i, let V^*(7j) be the set of nodes of the j-th component of 
C/i and let (7^ ) the set of nodes of the j-th component 
of C/i. The parameter f{gi) is evaluated averaging the 
individual contributions of the different connected com- 
ponents of gi or g.i according to the number of nodes 
they have against \Vi\ or \Vi\, the set of all nodes of CJ^ 
and gi, respectively, leading to: 

/(^^)^^ E l^'(7fe)l/(7fe) and 



/(^^) ^ ^ E |l^X7.)l/(7fc)- 



r(e. 



We impose, for both mathematical and conceptual con- 
sistency, that: 

{max{HUgc),Hf{gc)} = 0) ^ {T{g) = 0). 

Furthermore, if Ec = 0, (i.e., the case where the graph 
consists of a single node): 



T{g) = 0. 



4. The coordinates of hierarchy for directed graphs 

These three coordinates enabled us to define a mor- 
phospace, O, in which every net will be represented by a 
point u{g) = {T{g),F{g),0{g)) in this space. The co- 
ordinates defining such space are thought of to grasp the 
essentials of hierarchy, taking as the canonical hierarchi- 
cal structure the one satisfying the postulates presented 
in section [Bj Additionally, as we have seen, comparison 
between coordinates also provides us valuable informa- 
tion concerning the whole structural organization of the 
net. For example, the relation between O and F enables 
us to estimate how many cyclic regions are in the graph 
and at what point of the causal flow they are located. In 
addition, they provide an estimator of relevance of coor- 
dinate T, since low values of O, F indicate us that the net 
is mainly cyclical and that the fecd-forward-like analysis 
provided by T is less relevant than in the case where O, F 
are close to 1. To test the consistency of this coordinate 
system, we will rigorously check the behavior in several 
extreme and paradigmatic cases: 

Lemma: G displays T{g) = F{g) = 0{g) = 1 if and 
only if g is a directed tree in which all the leafs have the 
same length, as defined in section \C 1\ 

Proof 

(=>) By deflnition, this graph is a DAG, therefore, 
F{g) — 0{g) = 1. Furthermore, all its nodes but the 
one in M have fci„ — 1, which means that 

(vg' e , gLiGc)-!} Ui^i' GL(gc)-i})iHb{g') - 0). 

(DIO) 

And, since all but the set of minimal nodes will display 
kout > 2, then: 

{yg' e {gi, , 5^00-1} \J^5i, > o), 

(Dll) 

and, thus. 



18 



1. Obtaining the node-weighted condensed graph 



Original 
graph 




condensation 




Node-weighted 

condensed 

graph 



Forward layering 


Backward layering 






• /=1 


/=1 



Treenness 

Computation of entropies by layers 



li{\og 1 + li\og 2+li\og 2) ^^(i^log 4 +5^log 4 +i^log 2), 




Hf(G)= 0.5 




Gi 



]i{\og 1) i/l(log 1) i^dog 1 + log 1) M( V-^\og 2+ i^log 2 
1 1 



Hf(Gi) = Hi,(Gi) = Hf(Gi) = 
Trenness computation 



Ht(Gi) = l 



7=7 



^(f(G) +S(f(G,)+f(G,))) 



l"f'"bf 7= l^(-0.66 + - l)=-0.55 

Feedforwardness Orderability 

® ® ® 




F=^( 3/5 +3/5 + 2/2)+ ^^(2/5 + 2/5)=0.56 



000000( 



0=— =0.43 



FIG. 16: The definition of the hierarchy- based morphospace Q,. Visual representation of the computation of the three axes of 
the morphospace: treeness, feedforwardness and orderability {TFO coordinates, u). 



LiGc) 



1; ^i(ec)-l 



}){u,s.^{Hf{g'),H,{g') = Hf{g')). 



r 



The above claims are also true for the whole graph g 
Therefore, all the 2L{gc) — 1 quotients of the type shown 
in equation ( D8 1 involved in the computation of T{g) will 
have a value value equal to 1, which leads to T{g) = 1. 

(<^) If F{g) = 0{g) = 1, then g is a DAG. To see that 
T{g) — 1, it is enough to realize that the graph must 
satisfy equations (DIO Dll ) and that this can only hap- 



pen if C/ is a directed tree in which all the leafs have the 
same length. We observe that, if the laves have not the 
same length, then a will emerge in some computation, 



lowering the value to T{g) < 1. 

Corollary: g displays T{g) = -l,F{g) = 0{g) = 1 
if and only if g is a directed tree in which all the leafs 
have the same length, as defined in section \C 1\ but with 
changing the direction of all arrows. 

This latter case would belong to the perfect antihier- 
achical system. A graph consisting in a big cycle will 
be the paradigmatic example of non-hierarchical system, 
showing thus T{g) = F{g) = 0{g) = 0. Finally, 
a graph consisting of a chain of nodes, or a DAG in 
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FIG. 17; Illustration of the kind of graphs living in the rnorphospace fi. We located different kinds of toy graphs in their 
corresponding TFO coordinates, u, to grasp the kind of structures one might expect to observe in the different regions of the 
rnorphospace. Numbers in red indicate u for every graph. 



which all the layers have the same size and all nodes 
the same connectivity will be the paradigmatic example 
of non-hierarchical but ordered system, showing T{Q) — 
0, F{g) = 0{g) = 1. An illustration of the gallery of 
conformation on the region of possible networks in the 
space is provided in figure (17). Notice that some re- 



gions of the morphospace cannot be occupied by the very 
definition of the measures, while others represent strange 
configurations. It is worth to note that it is mathemat- 
ically possible to apply this formalism to networks con- 
taining more than one connected component. However, 
meaningful information can only be obtained from the 
study of a single component since it represents a unique 
causal structure. 



In figure ( 16 ) we detail the computation of the hier- 
archy coordinates of a given graph. From the original 
directed graph, the process of strongly connected compo- 
nent (see) detection and condensation give rise a node- 
weighted condensed graph by from the original graph. 
From this graph, the three coordinates are computed. 

All the theoretical work is thus finished. From now 
on, we will apply such machinery to the analysis of both 
model and real networks. 



Appendix E: Analysis of Networks 

In this section we systematically study the location of 
both real and model networks in the morphospace defined 
by the three coordinates T{g),F{g),0{g). We begin 
by studying classical models of random graphs, namely, 
the Erdos-Rcnyi (ER) |4], the Barabasi- Albert prefer- 
ential attachment (BA) [S] and the Callaway-Hopcroft- 
Kleinberg-Newman-Strogatz uniform attachment (Call- 



away) [S]. Then, we evaluate the hierarchy coordinates 
of 125 real networks belonging to different systems. The 
analysis of real nets is completed by confronting the ob- 
tained results with the exploration of the hierarchy co- 
ordinates of their corresponding randomized ensembles. 
We use two randomization methods which are exposed in 
detail. Finally we present an in silico experiment of evo- 
lution of networks inside the morphospace. Such analysis 
enables us to explore the accessibility of the space of pos- 
sible configurations and sheds light on what is likely to 
observe in real systems. 



1. Model Networks 

Analysis of model networks obtained using the directed 
versions of ER, BA and Callaway models show similar 
results. In all cases, model networks having high connec- 
tivities aggregate around the region defined by the rect- 
angle T = (-1, 1), F « (0, 0.5) O « 0. The high variance 
on the value of T can be explained by the fact that C/c 
is very small, since high connectivity produces a drastic 
process of condensation. It is easy to see that a network 
with a handful of nodes is more likely to give extreme 
values than large networks. Departing from the value 
O displayed at high connectivities, as connectivity 
decreases, O starts to increase but nets still aggregate 
around T « and F ^ 0. For low connectivities, T val- 
ues tend to be closer to T = but both F and O change 
their behavior, increasing their values to reach the region 
around the point T = Q,F = 1,0 = 1. This latter situ- 
ation can be explained by the low presence -or complete 
absence- of cycles in networks having small connectivi- 
ties. Actually, the region T = 0, F = 1, O = 1 is only 
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occupied by DAGs. Results are shown in figure (2b-c) 
of the main text Hierarchy in complex systems: the 
possible and the actual. 

Below we proceed to detail the exact parameters of the 
numerical experiments. 

Directed ER graphs.- We begin with the nets gen- 
erated by the directed version of the Erdos Renyi model. 
Such a model consists of two parameters, \V\ and p, that 
specify the size of the network {\V\) and the probabil- 
ity of connecting any pair by an arc in any direction 
(p). For \V\ ^ 0, this model does not ensure connect- 
edness until values of p ~ log(|F|)/|y| [7]. An ensemble 
of 35, 760 graphs with three sizes {\V\ = {100, 250, 500}) 
was created following the following procedure: We create 
64 replicas for every graph size starting from Pmax =0.1. 
We repeated this process reducing p in a step size of 
0.0005 until a value of p able to produce graphs with 
a fraction of vertices belonging to the giant connected 
component larger than 97.5%. The resulting ensemble 
encompasses a range of 100 < (fc) < 2.51. 

Directed Callaway graphs.- The uniform attach- 
ment [6^ is an iterative model consisting of three param- 
eters, namely mg, m and i. Starting on a fully connected 
set of mo nodes, during i iterations a new node is added 
to the network. Such node is linked to m existing nodes. 
These nodes can be selected with equal probability, i.e. 
there is a uniform attachment. Once the iterative process 
is finished, a random uniform orientation that converts 
links into directed arcs is performed in order to obtain a 
directed graph. An ensemble of 768 directed graphs was 
generated considering 64 replicas of networks for three 
different sizes(|y| = {100,250,500}) and four averages 
degrees ((fc) = {2,4,6,8}). Average degrees correspond 
with m = {1, 2, 3, 4}; in each case, the seed is a clique of 
mo — ni nodes. 

Directed BA graphs.- The preferential attachment 
[5] is an iterative model consisting of three parameters, 
namely mo, m and i. Starting on a fully connected set of 
mo nodes, during i iterations a new node is added to the 
network. Such node is linked to m of the existing nodes. 
These nodes are subject to be selected with a probability 
proportional to their degree, i.e. there is a preferential 
attachment. Once the iterative process is finished, a ran- 
dom uniform orientation that converts links into directed 
arcs is performed in order to obtain a directed graph. An 
ensemble of 768 directed graphs containing 64 replicas of 
networks for three different sizes(|y| = {100,250,500}) 
and four averages degrees {{k) = {2,4,6,8}, correspond- 
ing to m = {1, 2, 3, 4}) was generated. As above, the seed 
is a clique of mo = m nodes in all cases. 



2. Real Networks 

After looking at the properties of standard models of 
random networks, we explored a collection of 125 net- 
works encompassing 13 types of systems obtained from 
real data. Figure Results are shown in figure (2d) of the 



main text Hierarchy in complex systems: the pos- 
sible and the actual, shows the resulting TFO mor- 
phospace for the real networks used in this work. To dis- 
cuss the relevance of the observed results, we confronted 
data of real networks against the one obtained from their 
randomized counterparts. Randomized ensembles were 
built using two different techniques of randomization. 
Below methods and results are detailed, after the net- 
work data set presentation. 



a. Network dataset 

C elegans cell lineage network. Label: Cellular 
in the Fig 2 of the article. Cell lineage network obtained 
from Worm database and pulished as analysed in 

C elegans neural network. Label: Neuronal in 
the Fig 2 of the article. A directed, weighted net- 
work representing the neural network of C. elegans. 
Data compiled by D. Watts and S. Strogatz and made 
available at the Mark Newman's website: http://www- 



lpersonal.umich.edu /mej n / netdat a/ used in [9 . Original 
experimental data taken from [10] . 

Metabolic networks. Label: Metabolisms in the 
Fig 2 of the article.. Metabolic network set is a selec- 
tion of 19 reaction-metabolite directed networks obtained 
from three different published papers for different organ- 
isms: 

Barabasi dat aset for E.coli, B. subtilis, S cerevisiae. 

http: / / www.nd.edu/networks /resources / metabolic / index.html 
from data used in |li) . 

Networks obtained from the database published in 
P^ . A selection of 5 multicellular animals: H. sapi- 
ens (hsa), M. musculus (mmu), D. melanogaster (dme), 
R. novergicus (rno) and C. elegans (eel). One multi- 
celular plant: A. thaliana (ath). Two unicellular funghi: 
S. cerevisiae (see), S. pombe (spo). And seven prokary- 
otes: P. Aeruginosa (pae), E. coli (eco), B subtilis (bsu). 
Mycoplasma genitalis (mge). Mycoplasma pneumoniae 
(mpn), Synechocystis sp. PCC6803 (syn) and Salmonella 
typhimurium (sty). 

The Edinburgh human metabolic network taken from 

m- 

Food Webs. Label: Food webs in the Fig 2 of the 
article. These networks describes the exchange flow from 
the donor to recipient compartments where not neces- 
sarily nodes represent single species. The food web set 
contains a set of 22 directed weighted graphs in their 
original format. In this article were taken as unweighted 
directed graphs. Food-webs were originally selected from 
the R.E. Ulanowicz's Collection from the Ecosystem Net- 
work Analysis site and from ATLSS - Network Analysis of 
Trophic Dynamics in South Florida Ecosystems and com- 
piled in the Pajek dataset. Source: http://vlado.fmf.uni-| 
Ij .si / pub / networks /data/bio / foodweb /foodweb.htm. 
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Gene regulatory networks (GRNs) and one ki- 
nase network. Label: GRNs and Kinase network re- 
spectively, in the Fig 2 of the article. GRNs set contains 
the networks from two different articles. The data con- 
sists of directed networks where nodes are genes and arcs 
between two genes captures the interaction of the respec- 
tive gene product -represented by the source node- over 
the regulatory region of the target gene -target node. The 
only one available kinase network at the present is rep- 
resented by a directed network where nodes are proteins 
and arcs represents the relation of phosphorylation of a 
protein -node source- over the a target protein -node tar- 
get. GRNs of S cerevisiae, E. coli, Mus musculus, Ratus 
novergicus, Homo sapiens and Micoplasma tuberculosis 
and one kinase network of S cerevisiae were taken from 
http://info.gersteinlab.org/Hierarchy and published in 
'mi . GRNs of~ cerevisiae^ E. coli, and Bacillus sub- 
tilis taken from |15j . 

Electronic circuits. Label: Electronic circuits in 
the Fig 2 of the article. Directed networks of 50 elec- 
tronic wiring compiled from ISCAS'89 and ITC'99 sets. 
Data from electronic circuits used in (16j . 



Ownership:EVAownership. Label: Ownership in 
the Fig 2 of the article. EVA is a multidisci- 
plinary research project combining information extrac- 
tion, information visualization, and social network anal- 
ysis techniques to bring greater transparency to the 
public disclosure of inter-relationships between cor- 
porations. This data corresponds with an owner- 
ship network with 6,726 relationships among 8,343 
companies. |http: / / denali.berkeley.edu/e va/ Data ob- 
tained form pajek database: 'http: / /vlado.fmf.uni^ 
Ij. si/pub/networks/data/ econ/Eva/Eva.htm K. Norlen, 
G. Lucas, M. Gebbie, and J. Chuang. EVA: Extraction, 
Visualization and Analysis of the Telecommunications 
and Media Ownership Network. Proceedings of Inter- 
national Telecommunications Society 14th Biennial Con- 
ference (ITS2002), Seoul Korea, August 2002. 

Blogspol. Label: Social in the Fig 2 of the ar- 
ticle. A directed network of hyperlinks between we- 
blogs on US politics, recorded in [19]. Data obtained 
from pajek database. Original data was taken by pajek 
database from the Mark Newman's dataset |http: / / www- 
personal. umich.edu/ ~mejn/_ net data/ 



Word corpora. Label: Word corpora in the Fig 2 of 
the article. Four directed lexical graph from different texts 
according the description in [IT]: "Angle's Wren Xmas 
Tale" (Paul Auster), "Frankenstein or the modern prom- 
eteous" (Mary Shelley)," Ulisses" (James Joyce), "Moby 
dick" (Herman Melville), "Black cat" (Edgar Allan Poe) 
and a fragment of the New York Times column. 



Scientific citation network. Label: Citations in 
the Fig 2 of the article. A directed and feedforward graph 
of citations networks. Nodes are articles and every article 
points to the articles in which it is cited. Papers that cite 
S. Milgram's 1967 Psychology Tod ay paper or use Small 
World. Taken from pajek dataset: http://vlado.fmf.uni- 



|lj . si/ pub / networks/data / 



Software networks or dependence networks. 

Label: Software in the Fig 2 of the article. Ten directed 
networks capturing the map of file dependencies in 
different software. In these networks every node is a file. 
If a file A is called from another file B, then B receives 
an arrow from A. Although is generally avoided in pro- 
gramming, one file can content more than one function 
-argument- and therefore it can be called for more than 
one purposes. This feature would explain the appearance 
of cycles in some circumstances. Eight networks were 
taken from [18| . Two additional networks were obtained 
from the Internet: r package dependencies: taken from 
http://csgillespie.wordpress.com/2011/03/23/gTaphical- 
display-of-r-packages-dependencies/ Java packages taken 
from: http://gd2006.Org/contest/details.php#java, 



Hitech. Label: Social in the Fig 2 of the article. 
The case is a small hi-tech computer firm which sells, 
installs, and maintains computer systems. The network 
contains the friendship ties among the employees, which 
were gathered by means of the question: Who do you 
consider to be a personal friend? A friendship choice 
(arc) is only included in the network if both persons in- 
volved acknowledge it. Taken from http: / /vlad o.fmf .uni-, 
Ij.si/pub/networks/data/esna/hiTech.htm, published in 
to 

ModMath. Label: Social in the Fig 2 of the article. 
Taken from Pajek dataset (http://vlado.fmf.uni-lj.si/pub 



/networks/data/esna/modMath.htm. This network con- 
cerns the diffusion of a new mathematic method in the 
1950s. This innovation was instigated by top mathemati- 
cians and sponsored by the National Science Foundation 
of the USA as well as the U.S. Department of Educa- 
tion. R.O. Carlson, Adoption of Educational Innovations 
(Eugene: University of Oregon, Center for the Advanced 
Study of Educational Administration, 1965, p. 19). 



PhD advisors 

Genealogies in the Fig 

Data taken from: |http://vlado.fmf.uni-lj.si/ 



(genealogy). Label: 
2 of the article. 



pub/networks/data/esna/CSPhD.htm. The net- 

work contains the ties between Ph.D. students and their 
advisors in theoretical computer science; each arc points 
from a supervisor to a student. The partition contains 
the (estimated) year in which the Ph.D. was obtained. 
Original author: David Johnson; maintained by Ian 
Parberry. The SIGACT Theoretical Computer Science 
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Genealogy, Last Updated July 22, 1996. Data taken as 
appeared in [8] 



b. Randomization Methods 

Method A: Preserving the undirected degree 
sequence and the component structure. In this 
randomization technique, we generate the ensemble of 
randomized graphs by taking as topological invariants 
the undirected degree sequence and the component struc- 
ture. It is worth to note that the component structure is 
not preserved using a standard method of graph random- 
ization under the so-called configuration model approach 

m- 

Given an undirected graph Q{V,E), \V\ — n, its undi- 
rected degree sequence, au{G), is the sequence of n integer 
numbers in which the z-th number depicts the undirected 
degree of node Vi, k(vi) in a given labeling vi,...,Vn of 
the nodes of the graph: 

Let us suppose that such a graph has a given component 
structure T{Q) — 71,..., 7^. The randomization method 
applied here keeps invariant both au{Q) and Y{Q)- This is 
performed using a directed version of the local swap algo- 
rithm [HI [22] . Local swap diverges from standard random 
rewiring algorithms because it keeps invariant the com- 
ponent structure of the graph. The figure below shows 
how we generate the randomized elements of the ensem- 
ble. Let 5 be a directed graph. We choose three links 
such that they form a chain of length three in Q^, figure 
( 18 left), i.e., a "c" structure, no matter the direction of 
the arrows. Then, we cross the links at the extremes of 



the chain, i.e., generating a " k" structure, figure (18 cen- 
ter). Finally, we flip the senses of the arrows at random 
and, thus, links a and b have transformed into a' and b' . 
If a' or b' previously existed, we abort this randomizing 
event and we restart looking at random for another IZ 
structure within the graph. 





FIG. 18: A local swap algorithm keeping the undirected de- 
gree sequence invariant. 



For every real network we generated an ensemble of 
100 replicas obtained after applying the randomization 
algorithm until we performed A\E\ link switches or 20|£^| 
trials. The latter applies generally if the net belongs to 
an ensemble having a few members, which can happen 
due to several reasons, mainly, if the net is too dense 
or if the net is very sparse, or if it has a non-standard 
component structure. Numerical details of the random- 
ization process for each real network are described in the 
Appendix. 

Method B: Preserving the directed degree se- 
quence and the component structure. Now we 

generate the ensemble of randomized graphs by taking 
as topological invariants the directed degree sequence 
and the component structure. Given a directed graph 
Q{V,E), \V\ = n, its directed degree sequence, a{G), is 
the sequence of n pairs of integer numbers in which the 
first number of the i-th pair depicts the m-degree of node 
and the second one depicts the out-degree of node Vi, in 
a given labeling vi, ...,Vn of the nodes of the graph: 



<y{G) = {kin{vi), ^onf (^1)) ■■■7 (J^iniyi) •! ^out 

{Vi))..., {kin{ '^71 ) 7 ^OUt (Vn))- 



r 



Let us suppose that such a graph has a given component 
structure T{Q) = 71,..., 7^. The randomization method 
applied here keeps invariant both (j{G) and r(^). The 
conservation of cr(C/) makes this randomization method 
slightly more restrictive than the other presented above. 
Specifically, we look for C structures in G to perform 
a Local Swap, but, now, the sense of the arrows matter. 
Indeed, the kind of IZ structures over which we can apply 
a local swap keeping the directed degree sequence of the 
graph are: 

The kind of IZ structures over which we cannot apply 
a directed Local Swap operation, if we want to keep the 
directed degree sequence of the graph invariant are: 



Figure (21 ) shows how we generate the randomized el- 
ements of the ensemble. First, we select at random a 



structure of the kind described in figure (19), figure (21 
left). Then, we cross the links at the extremes of th( 
chain, i.e., generating a "k" structure (21 center). Fi- 
nally, we check if new links a' and 6' previously existed. 
If so, we abort this rewiring event and we restart an- 
other one looking at random for structures like the ones 
described in figure (19). 



For every real network we generated an ensemble of 
100 replicas obtained after applying the randomization 
algorithm until we performed 4|i?| link switches or 20|£^| 
trials, if the net belongs to an ensemble having a few 
members. 
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FIG. 19: Structures over which a directed Local Swap oper- 
ation can be apphed keeping the directed degree sequence of 
the graph invariant. 



FIG. 21: A local swap algorithm keeping the directed degree 
sequence invariant. 



i 



Randomization method a 



FIG. 20: Structures over which a directed Local Swap opera- 
tion cannot be applied in this case. 



c. Confronting real data with their randomized counterparts 




|1 u i' i ' 



Network randomization offers a quantitative picture of 
how far real networks are from their respective randomly 
generated ensembles whose correlations are eliminated by 
an iterative process of arc rewiring. Comparison is given 
by contrasting the TFO coordinates, u, of a particular 
real network with a distribution of values obtained from 
a randomized ensemble represented by its percentiles in 
a box plot fashion. This chart allows us to elude any 
assumption on the statistical behavior of the ensembles 
of randomized graphs. Figure (22) shows the compari- 



son between real and random ensemble data, the latter 
generated according methods a and b. Studied networks 
appear clustered by types, and every type of network has 
an associated color. Numbers, network labels, colors and 
other additional network quantifiers are detailed in the 
supplemental information II. In the figure, real network 
values are represented in charts by a grey square. The 
distribution of TFO values for every generated ensembles 
is represented in a box plot. Percentiles 25"^, and 75**^ 
are represented by an empty blue line box while the per- 
centile 50**^ is depicted with a red line inside. Whiskers 
show the values within percentiles lO"^ and 90'^. Finally, 
red crosses display the values out for the whisker range. 



In figure (22) looking at the position of real networks, 
they fall out from the percentile 50*^ of their randomized 
ensembles in the space. Therefore, we can argue that 
real networks of under scrutiny have not a representative 
graphical configuration for their degree sequence. This 
indicates that, although real networks tend to live in the 
null model regions as we saw earlier, system's constraints 
confer differences in the pattern of connections that dis- 
tinguish real networks from their randomized counter- 
parts. This difference is more accentuated for method a 
(the directed degree sequence is not conserved) than for 
method b, mainly due to the fact that the latter imposes 
a more dramatic restriction in the graphical configura- 
tion. 
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FIG. 22: Comparison of real networks with randomized en- 
sembles according methods a (upper part) and b (bottom 
part). In two groups of three charts, the figure displays three 
box plots for every randomization method showing the com- 
parison of T, F, O values of the real value with its respective 
randomized ensemble. X axis labels networks ordered by net- 
work ID. 



In relation to T coordinate, graphs tend to show T 
values close to zero. However, there are exceptions. Pos- 
itive biased values are observed in food webs and most of 
GRNs. Such nets are thus hierarchic (T > 0), in terms 
of how the causal flow is organized. On the contrary, 
electronic circuits show an anti hierarchical configuration 
(T < 0), with a tendency to occupy slightly negative, but 
significant, values of T. Looking at F and O values, real 
networks are generally far from the whiskers of the en- 
semble, when confronting real data against the ensemble 
obtained using the randomization method a -where di- 
rected degree sequence is preserved. In most cases real 
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data appears to avoid the general cyclic character ob- 
served in randomized networks. This trend seems to be 
reduced when the directed degree sequence is conserved 
(randomization method b). This points to the conclu- 
sion that not local correlation among connectivities but 
the pattern of inputs and outputs would explain part of 
TFO values. In general, such a trend presents a larger 
cyclic character in randomizations than in real networks. 
However, we can observe two exceptions to this general 
behavior. Looking at the O for randomized ensembles 
obtained through method a, metabolisms and word cor- 
pora show a higher cyclic character than their respective 
randomized ensembles. However, when compared to ran- 
domized ensembles obtained using method b, the trend 
is similar to the common values observed in randomized 
versions. As it happened in the previous cases, the range 
occupied by real and randomized data seems to be fairly 
justified in most cases by directed degree sequence. 



3. The accessibility of the Morphospace: Evolution 

Now we shall concern ourselves to an in silico ex- 
periment based on an evolutionary algorithm applied 
to graph population variation. So far we studied how 
both real nets and their corresponding random ensem- 
bles cover the different regions of the morphospace fl. 
The numerical experiment presented below will give us 
important information on how accessible are the different 
regions of the morphospace if we impose selective pres- 
sures related to specific values of hierarchy coordinates 
over the possible graphs. In raw words, we let evolve a 
population of graphs taking as selective pressure the dis- 
tance of such graphs to a given target point of fl. For 
a given graph, the closer it is to this target point, the 
higher is its selective value. At the global, evolutionary 
level, what we want to evaluate is how long it takes is to 
reach a given target point: The more it takes the pop- 
ulation to reach it, the more inaccessible is this region 
considered. Results are shown in figure (4) of the main 
text Hierarchy in complex systems: the possible 
and the actual. 

Let us detail the experiment. Using a method de- 
scribed in [23] we quantified the accessibility of a grid 
of points in the defined 3D morphospace 17. This exper- 
iment has a high computational cost, thus, we covered 
r2 by a grid containing 100 points in three cuts, corre- 
sponding to three TF planes at three different values 
of orderability. The three planes are the ones defined by 
O = {0.15, 0.5, 0.85}. For each decided O value, along the 
TF plane, minimum and maximum coordinates were dis- 
tanced a value of 0.05 distance units'33l from the bound- 
aries and all points of the grid were equally distanced at 
0.10. Each of the resulting TFO points was used as a 
target to be reached or approached by a euclidean dis- 
tance in the 3D morphospace. We say that to reach a 
given target TFO point , u*, is done when at least one 
graph of an evolving population acquires a TFO value 



u such that the euclidean distance between u and u* is 
smaller than 0.05 distance units. The initial condition of 
this algorithm is defined by a population of graphs P (al- 
ways single connected components), being |P| = 25, each 
one made of |y| — 25 nodes, generated following the di- 
rected Erdos Renyi model with p — 0.08 (see section E 1 
for model details). 

At each evolutionary iteration (generation) we perform 
the following steps: 

• Computation of the TFO values u{Q) for all Q £ P. 

• For each graph Q of P, we compute the euclidean 
distance (i(u((C/), u*) between vl{Q) and desired tar- 
get point u* e il. 

• We calculated the average of the euclidean distance 
over the population {d{u{Q),VL*)) . 

• We applied the selection criterion: we eliminate 

from the P those graphs satisfying d{u{{Q)^\i*) > 
{d{u{g),u*)). The set of survivals is P' = {G e P : 
d(u(5),u*) < (d(u(g),u*))}. 

• We create a sequence spi picking at random and 
copying elements from P' . We perform this op- 
eration |P\ P'\ times -therefore, the sequence has 
length £(sp/) = |P\P'|, and repetitions are allowed. 

• We apply the randomization operator over all 
graphs of the sequence, keeping the graphs of P' 
invariant. After this step, together with the un- 
changed graphs excluded from the rewiring process, 
we obtaining the new generation of P. 

This evolutionary algorithm is then performed a number 
of 1000 generations. The randomization operator is ap- 
plied to every graph of the sequence sp/ and consists in 
two steps: 

• Random addition from one up to three arcs. 

• Random deletion from one up to three arcs of a 
graph. This step must satisfy the condition of pre- 
serving the single connected component. If not, 
randomization event is aborted and new arcs are 
chosen to be removed. 

There is one situation defined by two conditions for 
which evolutionary algorithm can be stacked: Given a 
generation, i) if none of the graphs become close enough 
to the target point and ii) none or all graphs are be- 
low the fitness mean of the population. Then, to avoid 
the evolutionary algorithm freezes, the following rule is 
applied: 

• Half of the population is chosen at random to be 
eliminated and replaced by the other half, giving 
rise to P' . 

For each target point of the grid, 250 replicas of the 
evolutionary experiment were performed. Each popula- 
tion pursuing a target point was allowed to evolve up 
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to 1000 generations. The average number of generations 
needed to reach each target point was used as an esti- 
mator of the accessibihty of that part of il. Those cases 
were the accessibihty value is exactly 1000 indicate that 
none of the graphs in none of the 250 processes was able 
to reach that target point -see figure (4) of the main text 
for results. In these cases, we expect that more than 1000 
generations are needed. 

A note of caution must be added to interpret the 
strength of the obtained results. It is worth to note that 
parameters of the evolutionary algorithm may be deter- 
minant in shaping the observed fitness landscape, thereby 
accelerating or slowing the convergence process to the 
desired target point. The aim of this experiment was 
just to shed light on how accessible different regions of 



n are, considering selection and evolution driving forces. 
In this context, we numerically answered this question 
choosing a combination of parameters which we consider 
to be enough representative to provide us relevant infor- 
mation. And indeed, thanks to this approximation, we 
showed that not all regions are equally accessible but, in- 
stead, some regions are fairly improbable to be achieved. 
Among other factors, the size of the network may play a 
crucial role to achieve certain regions of the morphospace, 
since, for example it is easier for small networks to display 
extreme configurations. 

The described results create a nice, global picture of 
what is possible and what is actually observed within the 
zoo of complex networks. This objective is achieved from 
a rigorous definition of hierarchy. 
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[31] In the forthcoming lines, M = M{Qc), M = m(5c), unless 

the contrary is indicated. 
[32] Note that i) we ruled out the contributions of ^/L{gc) 
and GuGc) because they contain no links by definition, 
and, therefore, would impact to the computation of the 
hierarchy without justification, and ii) consistently, we 
compute the average between 2L{Qc) — 1 objects, the sie 



of the set W(C/) defined in equation (C6l. 
[33] Given two points m = (Ti, Fi, Oi), ua = (T2,J2,02) G 
Jl, the distance is computed using the standard euclidean 
norm, namely: 

d(ui,U2) = - T2Y + (Fi - F2)2 + (Oi - 02)2. 



